Implement implicit Euler and Crank-Nicolson methods.

Use `x` for the solution vector in each solver as temperature at the next time step, instead of `T`.

- forward Euler -> implicit Euler
- backward Euler -> explicit Euler

Previously, boundary conditions were applied after solving the tridiagonal system, which resulted in incorrect temperature values. This commit fixes the implementation by properly incorporating boundary conditions into the matrix coefficients before solving.

Changes:
  - Isothermal BC: Set matrix row to enforce T[i] = T_iso directly
  - Insulation BC: Use ghost point method to enforce ∂T/∂z = 0
  - Radiation BC: Keep as post-processing due to nonlinearity
  - Apply same fix to both Implicit Euler and Crank-Nicolson methods

  Add comprehensive boundary condition tests:
  - Test isothermal BC with non-zero temperatures (T_upper=100, T_lower=50)
  - Test insulation BC at lower boundary
  - Verify that boundary conditions are correctly satisfied
  - Add tests to main test suite in runtests.jl

This ensures that the implicit methods work correctly for all types of boundary conditions, not just the special case where T=0 at boundaries.

- Zero-conductivity case -> Calculate the radiative equilibrium temperature.
- Non-zero conductivity case -> Call a numerical solution method to calculate the temperature.

- Replace the Unicode middle dot (・) with ASCII caret (^) in the docstring of `heat_conduction_analytical.jl` to fix the "invalid escape sequence" error in Julia 1.6.
- This ensures compatibility with Julia 1.6 while maintaining the visual clarity of the ASCII art diagram.

…ction_analytical.jl`.